Due to the recent progress in experiments on trapped ultracold atoms,^[1–7] the nonequilibrium dynamics of isolated quantum systems has become a major subject of study in many-body physics.^[1–11] On one hand, real time evolutions of density profiles can be detected in experiments by the in situ imaging.^[12] On the other hand, advances in computational techniques allow us to calculate time evolutions for many-body systems and study the underlying physics with very high accuracy.^[13–20] Among those studies on nonequilibrium dynamics, there has been growing interest in studying the quantum transport phenomena in ultracold atoms. A mass current has been induced by a displacement or distortion of trap potential,^[12,21,22] a quench or imbalance of density or interaction,^[23–26] or changing the complex hopping for a system in an optical lattice with an artificial gauge field.^[27,28]

Pioneering experiments on generating artificial gauge fields in trapped ultracold atomic gases^[29–33] have paved the way to study some fundamental aspects of many-body physics. In particular, the emulation of synthetic electric and magnetic fields is crucial to extend their proven quantum simulation abilities further, e.g., to quantum Hall physics, topological insulators.^[34] In one dimension, when ultracold atoms are loaded in optical lattices, effective magnetic fields can be generated by combining a radio-frequency field and a Raman field,^[35] or by temporally modulating the optical lattice with designed patterns.^[36] Through Peierls substitution^[37] the hopping coefficient is generalized to a complex number , where amplitude t and phase can be tuned separately in experiments. Theoretically by ramping up the phase in a spatially uniform fashion, a finite quasi-steady-state current (QSSC) ensues from the nonequilibrium dynamics of noninteracting spinless fermions.^[27] The entanglement entropy does not change when the QSSC lasts.

In this paper, we will study the quench dynamics of noninteracting atoms loaded in one-dimensional (1D) optical lattices with artificial gauge fields. Initially we prepare the open finite-size system in the ground state. After suddenly changing the hopping coefficient (mainly the phase ) the system is driven out of equilibrium. Although the QSSC has been found in a very similar system,^[27] the nonequilibrium dynamics has not been studied in detail before. Here we study time evolutions of density and momentum distributions, and detailed properties of the mass current. The dynamics are compared for systems with different particle number, system size, statistics, and phase difference in hopping coefficient. Although the hard-wall confining potential, corresponding to the open boundary condition, has been recently realized in an ultracold atomic system,^[38] usually there is a harmonic trap. So the dynamics in trapped systems are also included. The dynamics can be experimentally tested with minimal modifications of available techniques.^[35,36]

This paper is organized as follows. In Section 2, we present the model and approach. In Section 3, we show time evolutions of density distribution, momentum distribution, and mass current at the center for finite uniform systems after suddenly changing the hopping coefficient. We also determine the relation between the amplitude of QSSC and parameters of the system both numerically and theoretically. In Section 4, we study the dynamics in trapped systems. A summary is presented in the last section.

We study the system of noninteracting ultracold atoms loaded in 1D optical lattices with artificial gauge fields, which are modeled by the following Hamiltonian 1 with an open boundary condition. ( ) is the creation (annihilation) operator of an atom at site j. L is the number of lattice sites. and are the amplitude and phase of the hopping coefficient, and they can be tuned separately. The system can be driven out of equilibrium by varying the hopping coefficient and here we study the sudden quench dynamics. Initially the system is in the ground state and modeled by the Hamiltonian with parameters and . After suddenly changing and to and , the system begins its quench dynamics. Notice that parameters and are site-independent.

In the nonequilibrium dynamics we will study time evolutions of density and momentum distributions. The density distribution is defined as 2 with the instantaneous state at time and the reduced single particle density matrix. The momentum distribution is defined as the Fourier transformation of the reduced single particle density matrix, In the dynamics, there is also a mass current flowing from the left half ( ) to the right half ( ) of the system. It is defined as with the particle number in the right half of the system. With the continuity equation, the mass current at the center can be calculated by 3 It can be measured experimentally by the method described in Refs. [12] and [24].

Since in a realistic ultracold atomic system there is usually an additional harmonic trap, we will also study the quench dynamics in trapped systems. The additional Hamiltonian for harmonic trap potential is 4 where V is the strength of harmonic trap and is the position of the trap center. The harmonic trap will remain unchanged in the dynamics.

We first study the dynamics of noninteracting spinless fermions in finite uniform systems. The total particle number is denoted by N and in numerical calculations the unit of time is set to 1. In Fig. 1(a) we show the contour plot of a typical time evolution of density distribution. The time evolution of the density distribution is independent of and , since all single particle states are independent of both of them. It depends on the phase difference , not on and themselves. This can be easily obtained from the diagonalized Hamiltonian and single particle wave functions.

Fig. 1.

	Figure Option ViewDownloadNew Window
	Fig. 1. (color online) Contour plots of time evolutions of density distribution (a) and momentum distribution (b) for a finite uniform system with 200 free fermions, L = 500, and .

For a system with different parameters the structure of the time evolution of the density distribution are the same. Initially the system is in the ground state, and the density distribution is uniform with small oscillations whose amplitude is proportional to due to the finite size effect. Here we have assumed that the filling factor does not satisfy or . After a sudden change of phase in the hopping coefficient, the density decreases at one end of the system and forms a shock wave, while at the other end the density increases and forms a rarefaction wave.^[28] As time goes on both density decreased and increased regions become larger and larger. Then these two parts meet together at some time and the density distribution becomes irregular after. Here we focus on the time interval before two parts meet together. In this time interval the density distribution consists of three plateaus: the plateau at the density decreased end (shock wave) with amplitude ; the plateau at the density increased end (rarefaction wave) with amplitude ; and the plateau at the center of the system with amplitude β. Fronts of shock and rarefaction waves move linearly with time; speeds of moving fronts only depend on the phase difference and the filling factor β. For most systems the absolute value of speeds is equal to 2. Particles move upwards (downwards) when the phase difference < 0 (> 0). Here has been assumed.

In Fig. 1(b) we show the time evolution of the momentum distribution for the same system as shown in Fig. 1(a). The time evolution of the momentum distribution only depends on , , and the filling factor. In the thermodynamic limit the spectrum of the system is . Then at time zero the momentum distribution is basically a typical free fermion momentum distribution centering around . After a sudden change of phase , except for the scattering caused by boundaries there is no other scattering process that can change momenta of particles. The momentum transfer for a particle scattered by boundary equals . Then at a finite time in the momentum distribution there is another block centering around . As time goes on the amplitude around decreases, while it increases around . The width of each block equals .

Next we study time evolutions of the noninteracting bosonic system, where we only have to study the single particle problem. Again time evolutions of density and momentum distributions are independent of and . A typical time evolution of the density distribution for a single particle system is shown in Fig. 2(a). Initially the system is in the ground state and the density distribution It is not uniform. After a sudden change of phase , the whole density distribution moves linearly (See Fig. 2(a)). However, a small part of the particles is scattered by boundaries and forms the rarefaction wave. Then the rarefaction wave interferes with the unscattered part and stripes appear. Before the interference becomes obvious, the linearly moving speed of the density distribution only depends on the phase difference . In order to determine this speed, we trace the position of maximum in the density distribution for a short period of time. After linearly fitting the data for the system with different phase differences, the speed is equal to .

Fig. 2.

	Figure Option ViewDownloadNew Window
	Fig. 2. (color online) Contour plots of time evolutions of the density distribution (a) and the momentum distribution (b) for a finite uniform system with one boson, L = 500, and . The value shown in the time evolution of the momentum distribution is .

In Fig. 2(b) we show the time evolution of the momentum distribution for the same single particle system as shown in Fig. 2(a). Initially the momentum distribution is basically a function centering around . Because of the scattering by boundaries, which is the only way to change particles' momenta, there will be another peak centering around at a finite time. However, for the bosonic system the density is much lower at two ends than at the center of the system. Then during the short time interval the proportion of particles scattered by boundaries is very small and so is the proportion of particles around . In order to show the scattered part clearly, we plot in Fig. 2(b).

Theoretically, for the dynamics of one particle governed by Hamiltonian (1), the single particle wavefunction, which can be expanded as , obeys 5 in the thermodynamic limit. Setting as the unit of time and using one can obtain the dispersion relation between wave vector (momentum) k and frequency (energy) : 6 The corresponding group velocity reads 7 For a single particle system, initially the momentum of the particle in ground state is in the thermodynamic limit. Then the group velocity of the particle , which is also the moving speed of the whole density distribution. For the noninteracting spinless fermionic system, the time evolution of a particle is not affected by others. Then the group velocities of these fermions are , . The maximal absolute group velocity of the fermionic system is 8 which is also the absolute moving speed of the shock and rarefaction wave fronts in the time evolution of the density distribution.

Next we study the time evolution of the mass current flowing from the left half to the right half of the system. Again we first study the noninteracting fermionic system. Time evolutions of currents for different finite uniform fermionic systems are shown in Fig. 3(a). The time evolution is independent of and , after setting as the unit of current. Finite quasi-steady-state currents (QSSCs) with constant amplitudes do exist (See also in Ref. [27]). Besides, there are oscillations on QSSCs whose amplitudes are proportional to . The QSSC ends when shock and rarefaction waves in the density distribution meet together (See also Fig. 1(a)). Then the duration of QSSC is linearly proportional to the system size and inversely proportional to the maximal absolute group velocity. The amplitude of QSSC depends on the phase difference and the filling factor β, not on , , N, and L themselves.

Fig. 3.

Figure Option ViewDownloadNew Window

Fig. 3. (color online) (a) Time evolutions of currents for finite uniform fermionic systems. for all systems. (b) Time evolutions of currents for systems with different ways of changing to . Corresponding are shown in the inset. Systems are with L = 500 and 200 fermions. (c) Time evolutions of currents for systems with a few particles. Systems are with L = 500 and . (d) The amplitude of QSSC versus the phase difference for fermionic systems with different filling factors. Inset in panel (d): versus the filling factor β for fermionic systems with L = 500 and .

In Fig. 3(b) we show time evolutions of currents for systems with different ways of changing to . The amplitude of QSSC is insensitive to the way of changing, while the duration of QSSC slightly depends on it. With a preliminary analysis of numerical data, it is easy to find out that the amplitude of QSSC decreases as the filling factor β moves away from the half filling 0.5. Then when the filling factor or , the amplitude of QSSC should be finite and small. But the QSSC does not exist in these cases because the amplitude of QSSC is comparable with the amplitude of oscillation caused by the finite size effect. Time evolutions of currents for fermionic systems with a few particles are shown in Fig. 3(c) and the QSSC indeed does not exist. To make the QSSC reappear, one can let , or increase the number of lattice sites. A typical time evolution of current for a single particle system is also shown in Fig. 3(c). No plateau can be found in the time evolution, namely, initially condensed noninteracting bosons do not support QSSC.

From the above discussions we know that only the fermionic system can support the QSSC. Now we are going to determine the relation between the amplitude of QSSC and system parameters. We first define an average of the mass current 9 where is the duration of QSSC. The average current can be numerically seen as the amplitude of QSSC. We already knew that the amplitude of QSSC only depends on the phase difference and the filling factor β. In Fig. 3(d) we show versus for systems with different filling factors, and in the inset we show versus β for systems with a given phase difference. Everyone can come up with the fact that these are sine functions. After fitting lots of data with sine functions, the amplitude of QSSC with .

Theoretically, under the open boundary condition the Hamiltonian (1) can be diagonalized into . The spectrum , and the corresponding single particle states Then the initial ground state can be written as In the dynamics, right after suddenly changing the hopping coefficient the mass current I at time reads 10 Here is the ground state of N−1 fermions. In the thermodynamic limit the Hamiltonian and the mass current operator commute. Then the mass current is conserved and is the amplitude of QSSC in the thermodynamic limit. Replacing summation with integration and letting , 11 in the thermodynamic limit.

For a fermionic system, we indeed have . With , and agree with each other very well for finite filling factors. The uncertainty of is proportional to because of the finite size effect (see Eq. (10)), the same behavior is for the amplitude of oscillations on the QSSC. For a fermionic system with filling factor (almost empty) or (almost full), according to Eq. (11). However, it is on the same order of the finite size effect, and no QSSC exists in the time evolution of the mass current. For a noninteracting bosonic system, the effective filling factor . Then as well as the finite size effect. No QSSC exists in the bosonic system.

In order to hold the atom cloud, there usually is a harmonic trap in a realistic ultracold system. For a trapped system the density distribution is not uniform. There can be a Mott insulator phase at the center, surrounded by superfluid phases in a trapped fermionic system. This inhomogeneity will cause different dynamics after a sudden change of hopping coefficient.

We first study the fermionic system. A typical time evolution of the density distribution for a trapped fermionic system is shown in Fig. 4(a). Here and later we set for the sake of simplicity. The density distribution is basically in Gaussian shape at time zero. And there can be a Mott insulator phase at the trap center for a system with sufficient fermions. Both of them have the same nonequilibrium behavior. After suddenly changing the phase , the whole density distribution moves linearly in a short time interval, with the direction and speed decided by the sign and value of the phase difference . Besides, particles have different group velocities at different positions, which depend on the local density and the phase difference (shown later). As time goes on the effect of the harmonic trap becomes obvious: the trap prevents particles from escaping and the kinetic energy of the system changes into the potential energy. At some time the moving direction of the density distribution changes and an oscillation begins. In Fig. 4(c) we show the long time evolution of the density distribution for the same system as shown in Fig. 4(a). A periodic oscillation shows up, and it is gradually damped by the lattice effect.^[39] The period of oscillation is determined only by the strength of the harmonic trap, and it is shorter for a tighter trap. The amplitude of oscillation is decided by the strength of the harmonic trap and also the phase difference.

Fig. 4.

Figure Option ViewDownloadNew Window

Fig. 4. (color online) Contour plots of time evolutions of the density distribution (a) and the momentum distribution (b) for a trapped system with 50 fermions, L = 600, , and . (c) The contour plot of the long time evolution of the density distribution for the same trapped system shown in panel (a). (d) The group velocity versus the local density for trapped fermionic systems with . Inset in panel (d): the group velocity versus the phase difference for trapped fermionic systems with 50 fermions, L = 600, and .

Next we study the dependence of the group velocity on the local density and the phase difference. Instead of studying group velocities at different positions, we study the group velocity at the position with maximal density for systems with different particle numbers. We trace the position of the maximal density in the time evolution of the density distribution for a short time interval. The maximal density basically does not change in this short time interval and we denote it by . The position changes linearly with time and the slope gives the group velocity. In Fig. 4(d) we show group velocities (dot) for systems with different and a fixed phase difference, and in the inset we show group velocities (dot) for systems with different phase differences and a fixed . After fitting the data for systems with different and phase differences, the group velocity reads 12 Results from this equation for the corresponding systems are also shown in Fig. 4(d) and the inset (red line). The linearly moving speed of the whole density distribution equals in the short time interval.

In Fig. 4(b) we show the time evolution of the momentum distribution for the same trapped fermionic system as shown in Fig. 4(a). The momentum distribution is also in Gaussian shape at time zero, centering around . After suddenly changing the phase , momenta of particles can change gradually due to the presence of the harmonic trap, which is different from the finite uniform system. The momentum distribution moves gradually with the direction, decided by the sign of the phase difference . At some time the momentum distribution reaches its extreme, where it centers around . Then a damped periodic oscillation begins, with the same period as in the time evolution of the density distribution. The amplitude of the oscillation is decided by the phase difference.

In the finite uniform system, the dynamics for a few particles and finite filling factor are very different. But they are very similar in the trapped systems. In Figs. 5(a) and 5(b) we show time evolutions of density and momentum distributions for a trapped system with only one particle. Here we only study the system with very weak harmonic trap and the local density is much smaller than 1. In the time evolution of the density distribution, the whole system basically has the same group velocity in a short time interval. The group velocity , which is the low density limit of Eq. (12).

Fig. 5.

	Figure Option ViewDownloadNew Window
	Fig. 5. (color online) Contour plots of time evolutions of the density distribution (a) and the momentum distribution (b) for a trapped system with one boson, L = 600, , and . (c) The current at the trap center versus time for a fermionic system with 50 fermions, L = 600, , and . Inset in panel (c): the local density at the trap center for the same system.

Finally, we study the mass current at the trap center. In a finite unform fermionic system, the constant QSSC ceases to exist right after the local density at the center changes. In the trapped system the density distribution is not uniform and the local density at the trap center usually changes with time (see the inset in Fig. 5(c)). Then usually there is no constant QSSC in a trapped system; however, the concept of the QSSC can be extended. In Fig. 5(c) we show a typical time evolution of mass current at the trap center (solid line). The amplitude of the mass current is indeed not a constant. In the uniform fermionic system the relation between the amplitude of the QSSC and the total density satisfies Eq. (11). Here in the trapped system, due to the inhomogeneity we thought it is reasonable to replace the total density with the time-dependent local density in Eq. (11). Then the equation depends on time and the QSSC at the trap center 13 Here is the time-dependent local density at the trap center. The corresponding and are shown in Fig. 5(c) (red dot) and in the inset. In a short time interval these two currents agree with each other very well. Then the QSSC can exist in the trapped fermionic system. The time interval for these two currents being much the same is mainly decided by the strength of the harmonic trap. According to Eqs. (12) and (13), one can obtain a QSSC at the trap center with constant amplitude, by tuning the local density at the trap center to 0.5 at time zero ( ).

We have studied the quench dynamics of noninteracting atoms loaded in 1D optical lattices with artificial gauge fields. In the presence of an artificial gauge field, the hopping coefficient can be a complex number . We assumed that initially the system is in the ground state with parameters and , and studied time evolutions of the density distribution, the momentum distribution, and the mass current at the center after suddenly changing the hopping to and . For a finite uniform fermionic system with finite filling factor β, the density distribution consists of three plateaus in a short time interval: the shock wave plateau with amplitude , the rarefaction wave plateau with amplitude at two ends, and the center plateau with amplitude β. The shock and rarefaction wave fronts move linearly toward the center with the speed equaling to the maximal absolute group velocity.

For a single boson system, the whole density distribution moves linearly in a short time interval, with the speed equaling to the group velocity . Both bosonic and fermionic systems have the similar time evolution of the momentum distribution. Initially the momentum distribution is a typical free fermion distribution centering around . Due to the scattering caused by boundaries, which is the only way to change momenta of particles, there is another block centering around at a finite time. In the time evolution of the mass current, only a fermionic system with finite filling factor can support a QSSC. The constant amplitude of QSSC . The lifetime of QSSC is determined by the system size and the maximal absolute group velocity. There is no QSSC in a bosonic system, because of the finite size effect.

In a trapped system the density distribution is not uniform. Time evolutions of density or momentum distributions for both fermionic and bosonic systems are similar. After suddenly changing the phase , the whole density distribution moves linearly in a short time interval with the speed equaling to . However, in a fermionic system particles at different positions have different group velocities . Long time evolutions of density and momentum distributions are damped oscillations.

What we studied here can be tested in the ultracold atomic system with minimal modifications of available experimental techniques.^[35,36] It is a very interesting and fundamental example of the transport phenomena and the nonequilibrium dynamics in ultracold atoms. Later we will study temperature effects and the system with interactions, which are much more complicated.^[28]